1 Traitement de données spatiales

1.1 Préparer son environement de travail

  • Créez un projet RStudio.
    Dans RStudio créez un projet et nommé le exercice. Il s’agit d’une bonne pratique qui vous facilitera la tâche. Cela améliore l’organisation et la portabilité de votre travail.
  • Dans le projet, créez un dossier data.
  • Téléchargez les données :
    • Un fichier contenant les polygones des communes d’Occitanie au format gpkg est disponible ici : Occitanie.gpkg. Il faut télécharger le fichier puis le placer dans le dossier data que vous venez de créer.
    • Un fichier contenant des données socio-démographiques sur les communes de France est disponible ici : base_cc_comparateur.xls. Il faut le télécharger et le mettre également dans le même dossier data.
  • Créez un fichier exo1.R et enregistrez le à la racine du projet.

Votre projet RStudio doit être structuré de la manière suivante :

├── data
│   └── base_cc_comparateur.xls
│   └── Occitanie.gpkg
├── exo1.R
└── exercice.Rproj

Vous devez enregistrer votre script régulièrement dans le fichier exo1.R.

Source du fichier Occitanie.gpkg : ADMIN EXPRESS COG édition 2019, IGN
Source du fichier base_cc_comparateur.xls : Base comparateur de territoires 2019, INSEE

1.2 Charger une couche géographique

  • Chargez la couche géographique des communes d’Occitanie en utilisant le package sf et la fonction st_read().
  • Vérifiez le système de projection avec st_crs(). S’agit-il de données projetée ? Si ce n’est pas le cas, transformez la couche des communes dans la projection française (Lambert 93, EPSG : 2154) avec la fonction st_transform().
library(sf)
## Linking to GEOS 3.7.1, GDAL 3.1.2, PROJ 7.1.0
occ_raw <- st_read(dsn = "data/Occitanie.gpkg", stringsAsFactors = FALSE)
## Reading layer `Occitanie' from data source `/home/tim/Documents/prz/carto_avec_r_exo/data/Occitanie.gpkg' using driver `GPKG'
## Simple feature collection with 4454 features and 11 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -0.3271723 ymin: 42.33292 xmax: 4.845565 ymax: 45.04557
## geographic CRS: WGS 84
st_crs(occ_raw)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World"],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
occ <- st_transform(x = occ_raw, crs = 2154)

1.3 Réaliser une sélection par attributs

  • Séléctionnez toutes les communes de la Haute-Garonne.
  • Enregistrez votre séléction dans un nouvel objet nommé com31
  • Afficher la nouvelle couche (en utilisant notament st_geometry())

Pour connaitre la liste de tous les noms ou code de région on peut utiliser la fonction unique().

unique(occ$INSEE_DEP)
##  [1] "32" "31" "82" "12" "81" "30" "11" "65" "46" "09" "34" "66" "48"
com31 <- occ[occ$INSEE_DEP == "31", ]
plot(st_geometry(com31))

1.4 Fusionner des entités

  • Fusionnez les communes de la région en un seul polygone (reg76), utilisez la fonction st_union ().
  • Créez la couche géograpique des départements de la région (dep76), utilisez la fonction aggregate() pour regrouper les polygones et calculer les sommes des populations communales.
  • Affichez les résultats (communes, départements et région)
reg76 <- st_union(occ)
dep76 <- aggregate(occ[,"POPULATION"], by = list(occ$INSEE_DEP), sum)
plot(st_geometry(occ), col = "lightblue1", border = "white", lwd = .5)
plot(st_geometry(dep76), col = NA, border = "lightblue2", lwd = 1, add = TRUE)
plot(reg76, col = NA, border = "lightblue3", lwd = 2, add = TRUE)

1.5 Créer une zone tampon

  • Créez une zone tampon d’une distance de 10 km autour des limites de la commune de Toulouse avec st_buffer(). Quel est le code INSEE de la commune de Toulouse?
toulouse <- com31[com31$INSEE_COM == "31555", ]
toulouse_buffer <- st_buffer(toulouse, 10000)
plot(st_geometry(com31), lwd = .5)
plot(st_geometry(toulouse), col = "darkblue", add = TRUE)
plot(st_geometry(toulouse_buffer), col = NA, lwd = 3, border = "red", 
     add = TRUE)

1.6 Réaliser une sélection par localisation

Déterminez quelles communes de la Haute-Garonne intersectent le buffer créé.

  • Utilisez la fonction st_intersects().
  • Inserez directement le resultat dans une variable nommée in_buffer dans l’objet com31.
com31$in_buffer <- st_intersects(x = com31, toulouse_buffer, sparse = FALSE)
head(com31, 2)
## Simple feature collection with 2 features and 12 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 504779 ymin: 6198834 xmax: 565151 ymax: 6270449
## projected CRS:  RGF93 / Lambert-93
##                         ID         STATUT INSEE_COM INSEE_ARR INSEE_DEP
## 4 COMMUNE_0000000009762857 Commune simple     31046         2        31
## 5 COMMUNE_0000000009760450 Commune simple     31547         1        31
##   INSEE_REG CODE_EPCI NOM_COM_M POPULATION TYPE NOM_COM
## 4        76 200072635     BAREN         12  COM   Baren
## 5        76 200068641   SEYSSES       8787  COM Seysses
##                             geom in_buffer
## 4 MULTIPOLYGON (((504930 6199...     FALSE
## 5 MULTIPOLYGON (((556904 6267...      TRUE

1.7 Afficher des couches géographiques

Affichez/superposez toutes les couches géographiques créees :

  • Les communes de la Haute-Garonne
  • Les communes intersectées par le buffer
  • Toulouse
  • Le buffer autour de toulouse
  • Les départements d’Occitanie
  • Les limites de l’Occitanie

Jouez sur les styles pour les différencier…

plot(reg76, col = NA, border = "grey50", lwd = 2, bg = "lightyellow")
plot(st_geometry(com31), col = "#aec8f2", border = "white", lwd = .5, 
     add = TRUE)
plot(st_geometry(com31[com31$in_buffer == TRUE, ]),col = "red", lwd = .5, 
     border = "white", add = TRUE)
plot(st_geometry(toulouse), col = "darkblue", border = "black", lwd = 1,
     add = TRUE)
plot(st_geometry(toulouse_buffer), col = NA, border = "black", lwd = 2, 
     lty = 2, add = TRUE)
plot(st_geometry(dep76), col = NA, border = "grey50", add = TRUE)

1.8 Créer une couche de points

Créer un objet sf (points) contenant la localisation de la préfécture de région à Toulouse.

  • Vous pouvez récupérer la longitude et la latitude sur Google Map ou OSM.
  • Créer un point avec st_points(), puis un objet sfc avec avec st_sfc() et finalement un objet sf avec st_sf().
# Sur OSM : https://www.openstreetmap.org/search?whereami=1&query=43.59825%2C1.45047#map=19/43.59825/1.45047
# Position de la préfecture 43.59825, 1.45047
library(sf)
pref_pt <- st_point(c(1.45047, 43.59825))
pref_sfc <- st_sfc(pref_pt, crs = (4326))
pref <- st_sf(name = "Préfecture", geometry = pref_sfc)
pref
## Simple feature collection with 1 feature and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 1.45047 ymin: 43.59825 xmax: 1.45047 ymax: 43.59825
## geographic CRS: WGS 84
##         name                 geometry
## 1 Préfecture POINT (1.45047 43.59825)

1.9 Calculer des distances

Calculez une matrice de distance entre la préfecture et les centroïdes des communes du département Haute-Garonne.
Ajouter cette distance dans une nouvelle colonne dist_pref dans com31. Pour cela, utiliser les fonction st_centroid() et st_distance().
N’oubliez pas de vérifier les projections utilisées avec st_crs(), au besoin modifiez les avec st_transform().

st_crs(pref)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World"],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
st_crs(com31)
## Coordinate Reference System:
##   User input: EPSG:2154 
##   wkt:
## PROJCRS["RGF93 / Lambert-93",
##     BASEGEOGCRS["RGF93",
##         DATUM["Reseau Geodesique Francais 1993",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4171]],
##     CONVERSION["Lambert-93",
##         METHOD["Lambert Conic Conformal (2SP)",
##             ID["EPSG",9802]],
##         PARAMETER["Latitude of false origin",46.5,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",3,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",49,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",44,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",700000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",6600000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["France"],
##         BBOX[41.15,-9.86,51.56,10.38]],
##     ID["EPSG",2154]]
identical(st_crs(pref), st_crs(com31))
## [1] FALSE
pref <- st_transform(pref, st_crs(com31))
identical(st_crs(pref), st_crs(com31))
## [1] TRUE
# Centroides des communes
com31_centro <- st_centroid(st_geometry(com31))
com31$dist_pref <- st_distance(com31_centro, pref)
head(com31, 2)
## Simple feature collection with 2 features and 13 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 504779 ymin: 6198834 xmax: 565151 ymax: 6270449
## projected CRS:  RGF93 / Lambert-93
##                         ID         STATUT INSEE_COM INSEE_ARR INSEE_DEP
## 4 COMMUNE_0000000009762857 Commune simple     31046         2        31
## 5 COMMUNE_0000000009760450 Commune simple     31547         1        31
##   INSEE_REG CODE_EPCI NOM_COM_M POPULATION TYPE NOM_COM
## 4        76 200072635     BAREN         12  COM   Baren
## 5        76 200068641   SEYSSES       8787  COM Seysses
##                             geom in_buffer     dist_pref
## 4 MULTIPOLYGON (((504930 6199...     FALSE 104829.70 [m]
## 5 MULTIPOLYGON (((556904 6267...      TRUE  17211.48 [m]

1.10 Cartographier une variable

Cartographiez distance de chaque commune du département à la préfecture avec la fonction plot().

# Le plus simple
plot(com31["dist_pref"])

# Avec un peu de paramétrage
plot(com31["dist_pref"], main = "Distance à la préfecture (en mètres)", 
     pal = hcl.colors(13,"Turku", rev = TRUE), border = NA,
     key.pos = 1,key.width = .15, key.length = .75, graticule = TRUE,
     reset = FALSE)
plot(st_geometry(pref), pch = 20, col = "red", add = TRUE )

1.11 Extraire des données OpenStreetMap

  • Sélectionnez les communes d’un autre département dans un nouvel objet nommé com+codeDuDépartement* (com46 par exemple).
  • Avec le package osmdata extrayez l’ensemble des restaurants présents dans le département.
  • Affichez les communes qui ne contiennent pas de restaurants.
  • Pour ces communes calculez la distance au restaurant le plus proche.
com46 <- occ[occ$INSEE_DEP == "46", ]
plot(st_geometry(com46))

library(osmdata)
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
# Définition d'une bounding box
q <- opq(bbox=st_bbox(st_transform(com46,4326)))
# Extraction des restaurants
res <- add_osm_feature(opq = q, key = 'amenity', value = "restaurant")
res.sf <- osmdata_sf(res)
res.sf.pts  <- res.sf$osm_points[!is.na(res.sf$osm_points$amenity),]
resto <- st_transform(res.sf.pts, st_crs(com46))
resto <- st_intersection(resto, st_geometry(com46))
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
# Affichage des restaurants
plot(st_geometry(com46), col="darkseagreen3", border="darkseagreen4")
plot(st_geometry(resto), add=TRUE, pch=20, col = "#330A5FFF", cex = 0.5)

# Compter les restaurants par communes
inter <- st_intersects(x = com46, y = resto)
com46$nresto <- sapply(inter, length)

# communes sans restos
com46noresto <- com46[com46$nresto==0, ]
plot(st_geometry(com46), col="darkseagreen3", border="darkseagreen4")
plot(st_geometry(com46noresto), col = 'red', add = TRUE)
plot(st_geometry(resto), add=TRUE, pch=20, col = "blue", cex = 1)

# index du restaurant le plus proche
index <- st_nearest_feature(x = st_centroid(com46noresto), 
                            y = resto)
## Warning in st_centroid.sf(com46noresto): st_centroid assumes attributes are
## constant over geometries of x
# distance au plus proche
com46noresto$dresto <- st_distance(x = st_centroid(com46noresto), 
                                   y = resto[index, ], 
                                   by_element = TRUE)
## Warning in st_centroid.sf(com46noresto): st_centroid assumes attributes are
## constant over geometries of x
# Affichage de la carte
plot(com46noresto['dresto'], reset = F, 
     main = "Distance au restaurant le plus proche")
plot(st_geometry(com46), col=NA, add= TRUE)
plot(st_geometry(resto), add=TRUE, pch=20, col = "red", cex = 1)


2 Cartographie thématique

2.1 Préparer son environement de travail

Dans le même projet RStudio créez un deuxième script nommé exo2.R.
Votre projet RStudio doit maintenant être structuré de la manière suivante :

├── data
│   └── base_cc_comparateur.xls
│   └── Occitanie.gpkg
├── exo1.R
├── exo2.R
└── exercice.Rproj

2.2 Charger des données statistiques

  • Charger la couche géographique des communes d’Occitanie (data/Occitanie.gpkg). Utilisez le package sf et la fonction st_read() pour importer les données.
  • Vérifiez le système de projection avec st_crs().
  • Charger les fichier de données data/base_cc_comparateur.xls fournie par l’INSEE. Utilisez le package readxl et la fonction read_excel() pour ouvrir la table de données correctement. Importer la table dans un objet nommé occ_df.
library(sf)
occ_raw <- st_read(dsn = "data/Occitanie.gpkg", stringsAsFactors = FALSE)
## Reading layer `Occitanie' from data source `/home/tim/Documents/prz/carto_avec_r_exo/data/Occitanie.gpkg' using driver `GPKG'
## Simple feature collection with 4454 features and 11 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -0.3271723 ymin: 42.33292 xmax: 4.845565 ymax: 45.04557
## geographic CRS: WGS 84
st_crs(occ_raw)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World"],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
occ <- st_transform(x = occ_raw, crs = 2154)
library(readxl)
occ_df <- read_excel(path = "data/base_cc_comparateur.xls", sheet = 1, skip = 5) 

head(occ_df, 2)
## # A tibble: 2 x 36
##   CODGEO LIBGEO REG   DEP   P16_POP P11_POP SUPERF NAIS1116 DECE1116 P16_MEN
##   <chr>  <chr>  <chr> <chr>   <dbl>   <dbl>  <dbl>    <dbl>    <dbl>   <dbl>
## 1 01001  L'Abe… 84    01        767     780  16.0        41       25     306
## 2 01002  L'Abe… 84    01        243     234   9.15       21        7     101
## # … with 26 more variables: NAISD18 <dbl>, DECESD18 <dbl>, P16_LOG <dbl>,
## #   P16_RP <dbl>, P16_RSECOCC <dbl>, P16_LOGVAC <dbl>, P16_RP_PROP <dbl>,
## #   NBMENFISC16 <dbl>, PIMP16 <dbl>, MED16 <dbl>, TP6016 <dbl>,
## #   P16_EMPLT <dbl>, P16_EMPLT_SAL <dbl>, P11_EMPLT <dbl>, P16_POP1564 <dbl>,
## #   P16_CHOM1564 <dbl>, P16_ACT1564 <dbl>, ETTOT15 <dbl>, ETAZ15 <dbl>,
## #   ETBE15 <dbl>, ETFZ15 <dbl>, ETGU15 <dbl>, ETGZ15 <dbl>, ETOQ15 <dbl>,
## #   ETTEF115 <dbl>, ETTEFP1015 <dbl>

2.3 Réaliser une sélection par attributs

Créez un nouvel objet sf à partir d’une séléction par attribut`:

  • Séléctionnez toutes les communes d’un seul département (ne choisissez pas la Haute-Garonne (31)).
  • Enregistrez votre séléction dans un nouvel objet nommé com+codeDuDépartement* (com46 par exemple)
  • Afficher la nouvelle couche (plot(st_geometry(...)).

Pour connaitre la liste de tous les noms ou code de région voir la fonction unique()

unique(occ$INSEE_DEP)
##  [1] "32" "31" "82" "12" "81" "30" "11" "65" "46" "09" "34" "66" "48"
com31 <- occ[occ$INSEE_DEP=="31", ]
plot(st_geometry(com31))

2.4 Réaliser une jointure

  • Joignez la table de données INSEE avec la couche gégographique des communes de votre département séléctionné. Utilisez la fonction merge(). Quel est l’identifiant commun entre la table géo et la table attributaire?

2.5 Réaliser des cartes thématiques

Réalisez huit cartes thématiques finalisées.

Chacune des cartes devra comprendre les éléments d’habillage et de mise en page nécessaires leur compréhension (titre, légende, sources…)

Parmi ces huits cartes deux échelles différentes (échelle de la région, échelle d’un département) et deux découpages territoriaux différents (les communes, les départements) doivent être utilisés.

Vous devez réaliser les huit types de cartes suivants :

  1. Une carte de localisation (avec carton de localisation si possible, des labels, legende…);
  2. Une carte représentant une variable quantitative absolue (un stock);
  3. Une carte représentant une variable quantitative relative (un ratio);
  4. Une carte combinant une variable quantitative absolue et une variable quantitative relative (un stock et un ratio);
  5. Une carte représentant une variable qualitive;
  6. Une carte utilisant une anamorphose. Il est possible de combiner les anamorphoses avec d’autres variables (qualitative ou quantitative relative);
  7. Une carte interactive des communes du département + un point sur la préfecture du département;
  8. Une carte utilisant une grille régulière ou représentant des discontinuités.

Exemple :

  • Une carte en symboles proportionels représentant la population des départements de la région (échelle régionale, découpage départementale, carte de stock.).
  • Une carte choroplète représentant l’évolution de la population des communes d’un département (échelle départementale, découpage communale, carte de ratio).

Nous vous conseillons d’utiliser le package cartography pour l’ensemble de vos réalisations. Néanmoins l’utilisation d’autres packages comme ggplot2 ou tmap est possible.

2.6 Évaluation

Vous devrez nous rendre pour le xx février 2020 un script permettant de construire ces huit cartes dans un fichier dont le nom est constitué de votre prénom et de votre nom (Prenom_Nom.R).
Votre script doit pouvoir fonctionner une fois inséré dans un projet RStudio structuré de la manière suivante :

├── data
│   ├── base_cc_comparateur.xls
│   └── Occitanie.gpkg
├── Prenom_Nom.R
└── exercice.Rproj

Vous serez évalués de la manière suivante :

  • 2 points par cartes (choix sémiologiques, traitement des données en amont…),
  • 4 points pour la qualité et la présentation du code (reproductibilité, commentaires, clareté).

2.7 Aides

2.7.1 Créer de nouvelles variables

Pour répondre à cet exercice vous devrez créer au moins deux nouvelles variables. Une variable quantitative relative :

Une variable qualitative (n’utilisez par cette variables dans votre exercice ) :

2.7.2 Créer un carton avec la position de votre département


  1. UMS RIATE, CNRS↩︎

  2. FR CIST, CNRS↩︎